********************************************************************************
*********** This program is written for the present value project **************
********************************************************************************
* 02. VAR Analysis Firm Level
********************************************************************************

 

********************************************************************************
* 02.00  Initialize STATA
********************************************************************************
clear all
set more off, perm

********************************************************************************
* 02.01  Settings
********************************************************************************
* vweight?  0 = No  1 = Yes
local vweight = 0
local vweight2 = 0

* Extra state variables in augmented VAR
local basic_vars = "r roe scale bm"
local extra_vars = "assetgrowth"  

* Sample period
local start = 1965																 
local end = 2022

* Subsample cutoff
local late = 1990

* Choose rho
local rho = 0.96 

* Size, MB and Vuolteenaho requirement (1: yes, 0: no)
local size_req = 1
local mb_req = 1
local vuol_req = 0 														

* Include error analysis in the run?
local erroranalysis = 1

* Include subsample VAR in the run?
local subsamples = 0

* Characteristics used in subsample analysis
if `subsamples'==1 {
	local sort_chars = "mb_L1 ME_L1 ROE_L1 Dur_MW_L1 year" 
	local subsample_header1 = "&\multicolumn{2}{c}{{\it Growth }}&\multicolumn{2}{c}{{\it Large }}&\multicolumn{2}{c}{{\it Profitable }}&\multicolumn{2}{c}{{\it High duration }}&\multicolumn{2}{c}{{\it Recent }}" 
	local subsample_header2 = "&\multicolumn{2}{c}{{\it Value }}&\multicolumn{2}{c}{{\it Small }}&\multicolumn{2}{c}{{\it Unprofitable }}&\multicolumn{2}{c}{{\it Low duration }}&\multicolumn{2}{c}{{\it Early }}" 
}


********************************************************************************
* 02.02  Apply data filters and save temporarily
********************************************************************************
* Load annual dataset
use "Data/Created/Annual Data 2021 with Bins", clear 

* Choose sample years
keep if year >= `start' & year <= `end' 

* Default requirements
keep if datareq==1 & ffdata_L1==1

* Optional requirements
if `size_req' == 1 {
	keep if size_L1 == 1
}
if `mb_req' == 1 {
	drop if MB_L1 < 1/100 | MB_L1 > 100 | MB < 1/100 | MB > 100
}	
if `vuol_req' == 1 {
	keep if vuolreq == 1
}  

* Divide samples into two
if `subsamples'==1 { 
	gen bin_year = 1
	replace bin_year = 10 if year >= `late'
}

* Save temporarily
tempfile data
save "`data'"

* Save permanently for the bootstrap code
*save "Data/created/Annual Data Feb2021 - for bootstrap", replace


********************************************************************************
* 02.03  Summary stats   (Table 02_03)
******************************************************************************** 

local count = 1
foreach var of varlist r roe scale mb {
	preserve
	
	* Equal-weight all years
	collapse (count) N = `var' (mean) mean = `var' (sd) sd = `var' (min) min = `var' ///
			 (p1) p1 = `var' (p25) p25 = `var' (p50) p50 = `var' (p75) p75 = `var' (p99) p99 = `var' (max) max = `var', by(year) 
	collapse (sum) N (mean) mean sd min p1 p25 p50 p75 p99 max
	
	* Save stats
	qui tostring N, force replace format("%15.0fc")
	qui tostring mean-max, force replace format("%5.3f")
	g var = "`var'"
	order var
	tempfile tmp`count'
	save "`tmp`count''"
	
	restore
	local count = `count' + 1
	}

* Append all stats 
u "`tmp1'", clear
forval i = 2/4 {
	append using "`tmp`i''"
	}
	
* Label variabels 
replace var = "$ r_{t}$" 		if var == "r"
replace var = "$ mb_{t}$" 		if var == "mb"
replace var = "$ roe_{t}$" 		if var == "roe"
replace var = "$ iva_{t}$" 		if var == "scale"  
replace var = "\addlinespace " + var if _n > 1

* Export to Latex 
listtab * using "Results/Tables/02_03_`end'.tex", ///								* COMMENTS: IMPOSE ALL REQUIREMENTS
		rstyle(tabular) replace ///															HERE AND MAKE SURE THE NUMBER OF 
		head("\begin{tabular}{l c c c c c c c c c c}" ///									OBSERVATIONS SAME ACROSS ROWS.
		"\midrule" ///
		"Variable & N & Mean & St. Dev. & Min & 1\% & 25\% & Median & 75\% & 99\% & Max \\" ///
		"\midrule") ///
		foot("\midrule" "\end{tabular}")	



********************************************************************************
* 02.04  Error comparison  (Table 02_04)
********************************************************************************
if `erroranalysis'==1 {

	use "`data'", clear

	* CKLP approximation
	gen e_CKLP = -mb_L1 -r + roe + scale + `rho'*mb
	gen compare_e_CKLP = .
	gen r_CKLP = -mb_L1 + roe + scale + `rho'*mb
	* Campbell-Shiller
	gen e_CS = -m_L1 -log(`rho')-(1-`rho')*log(1/`rho'-1) -r + (1-`rho')*dfirm + `rho'*m_adj
	gen compare_e_CS = cond(abs(e_CKLP) < abs(e_CS), 100, 0)
	gen r_CS = -m_L1 -log(`rho')-(1-`rho')*log(1/`rho'-1) + (1-`rho')*dfirm + `rho'*m_adj
	* Vuolteenaho
	gen e_V = -mb_L1 -r + roe + `rho'*mb 
	gen compare_e_V = cond(abs(e_CKLP) < abs(e_V), 100, 0)
	gen r_V = -mb_L1 + roe + `rho'*mb 
	 
	* AAPL and IBM figures
	foreach gvkey in "001690" "006066" {
		twoway (line r year if gvkey=="`gvkey'", lwidth(vthick) lcolor(gs10)) ///
			(line r_CKLP year if gvkey=="`gvkey'", lwidth(medthick) lpatter(-) lcolor(red)  ///
			  graphregion(color(white) lcolor(white)) ///
			   bgcolor(white) ytitle("Actual or approximated return", size(large)) ///
			   xtitle("",)  ylabel(, labsize(large)) ///
			   xlabel(, labsize(large)) legend(label(1 "Actual log return") size(large) ///
			   label(2 "Approximation")))
		gr export "Results/Figures/02_04_`gvkey'.pdf", replace as(pdf) 
	} 
	
	* Work on the table 
	local count = 1
	foreach var of varlist e_CKLP e_CS e_V {
		preserve
		
		* Equal-weight all years
		collapse (count) N = `var' (mean) compare = compare_`var' (mean) mean = `var' (sd) sd = `var' (min) min = `var' ///
					 (p1) p1 = `var' (p25) p25 = `var' (p50) p50 = `var' (p75) p75 = `var' (p99) p99 = `var' (max) max = `var', by(year) 
		collapse (sum) N (mean) compare mean sd min p1 p25 p50 p75 p99 max
		qui tostring N, force replace format("%15.0fc")
		qui tostring compare, force replace format("%15.1fc")
			replace compare = compare + "\%"
		qui tostring mean-max, force replace format("%5.3f")
		g var = "`var'"  
		order var
		tempfile tmp`count'
		save "`tmp`count''"
		
		* Error comparison for large cap
		restore
		preserve 
		collapse (mean) compare_large = compare_`var' if bin_ME_L1==10, by(year) 
		collapse compare_large
		qui tostring compare_large, force replace format("%15.1fc")
			replace compare_large = compare_large + "\%"
		g var = "`var'"
		merge 1:1 var using "`tmp`count''", nogen
		order var N mean-max compare compare_large
		tempfile tmp`count'
		save "`tmp`count''" 
		
		* Restore data and update count
		restore
		local count = `count' + 1
		}
		
	* Append all stats 
	preserve	
	u "`tmp1'", clear
	forval i = 2/3 {
		append using "`tmp`i''"
		}    
	replace compare = "" 			if var=="e_CKLP"
	replace compare_large = "" 		if var=="e_CKLP"
		
	* Label variabels 
	replace var = "CKLP" 					if var == "e_CKLP"
	replace var = "Campbell-Shiller" 		if var == "e_CS"
	replace var = "Vuolteenaho w/o CS adj" 			if var == "e_V"
	replace var = "\addlinespace " + var 	if _n > 1

	* To save space, drop N
	drop N

	* Export to Latex 
	listtab * using "Results/Tables/02_04_`end'.tex", ///								* COMMENTS: IMPOSE ALL REQUIREMENTS
			rstyle(tabular) replace ///															HERE AND MAKE SURE THE NUMBER OF 
			head("\begin{tabular}{l c c c c c c c c c c c}" ///									OBSERVATIONS SAME ACROSS ROWS.
			"\midrule" ///
			"& \multicolumn{9}{c}{Approximation error $ e_t $} & \multicolumn{2}{c}{$|e_t|>|e^{CKLP}_{t}|$} \\" ///
			"\cmidrule(lr){2-10} \cmidrule(lr){11-12}" ///
			"Model & Mean & St. Dev. & Min & 1\% & 25\% & Median & 75\% & 99\% & Max & All stocks & Top size decile \\" ///
			"\midrule") ///
			foot("\midrule" "\end{tabular}")
	restore		
	
}
	
	
********************************************************************************
* 02.05  Begin VAR analysis for spec==1   
********************************************************************************

local spec = 1 // specification: 1 (one lag, 4 variables), 2(one lag, include extra variables)  

	use "`data'", clear

	* Number of lags to use
	local lags = 1
 
	* Weight on each observation ((#firms in the year)^(-1) * possible value weight)
	tempfile tmp
	save "`tmp'" 
	collapse (count) firms = mb (sum) ME_L1, by(year)
	ren ME_L1 sumME_L1
	merge 1:m year using "`tmp'", nogen
	g weight = 1/firms
	if `vweight'==1 {
		replace weight = ME_L1/sumME_L1
	}
	gen valueweight = ME_L1/sumME_L1
 
	* Choose state variables	
	local state_vars = "`basic_vars'"  
		 
	* Choose all variables including lags   % CHANGED LOOP ORDER to vars inside lags, so regressors will be ordered like vars when lags>1
	local all_vars = "`state_vars'" 
	forval i = 1/`lags' {
		foreach var in `state_vars' {
			if "`var'"=="bm" & `i'<=1 {
				local all_vars = "`all_vars' `var'_L`i'" 
			}
			else if inlist("`var'","roe","scale","sgrowth","gamma") & `i'<=2 {
				local all_vars = "`all_vars' `var'_L`i'" 
			}
			else if inlist("`var'","r","assetgrowth") {
				local all_vars = "`all_vars' `var'_L`i'"  
			}
		}
	}	 
		
	* Count the total number of variables   
	local num_vars = 0 
	foreach var in `state_vars' {
		local num_vars = `num_vars' + `lags' 
	} 
	
	* Define matrix to store the results
	matrix B = J(1,`num_vars',.)													// Coefficient matrix
	matrix SE = J(1,`num_vars',.)													// Standard error matrix
	matrix Rsq = J(1,1,.)
 
	* Define regressors and set of variables
	local regressors = "" 														
	foreach var in `state_vars' {
		local regressors = "`regressors' `var'_L1" 
	}				
 
 
****************************************************************************
* 02.06  Regress state variables  (Table 02_06_`vweight')
**************************************************************************** 
	local i = 1
	foreach var of varlist `state_vars' {
		 
		* Regression
		reghdfe `var' `regressors' [aweight = weight], nocon absorb(year) cl(year firm) resid
			 
		* Compute error terms				
		predict u_`var', residual
					 
		* Store coefficients and standard errors						
		matrix B = (B \ (e(b)))
		matrix V = vecdiag(e(V))
		matmap V V, map(sqrt(@))
		matrix SE = (SE \ V)
		matrix Rsq = (Rsq \ (e(r2_within)))
	 
		local i = `i' + 1	
	} 
	 
	* Remove first row
	mat B = B[2...,1...]
	mat Rsq = Rsq[2...,1...]
	mat SE = SE[2...,1...]
	  
	* Save data on residuals
	tempfile resid
	save "`resid'"
	
	* Create a table of coefficients
	* Import stored coefficients and standard errors 
	clear 
	svmat B
	svmat Rsq
	svmat SE
	
	ren Rsq* Rsq
	
	* Standard erros  
	forval i = 1/`num_vars'{ 
		qui tostring SE`i', force replace format("%5.3f")
		qui replace SE`i' = "(" + SE`i' + ")"
	} 
	
	* Make a column for variable names	
	g col = ""
	local j = 1
	foreach var in `state_vars' {
		if "`var'"=="scale" {
			replace col = "$ 0`j'iva_{t} $" if [_n]==`j'
		}
		else {
			replace col = "$ 0`j'`var'_{t} $" if [_n]==`j'
		}
		local j = `j' + 1
	}  
	replace col = "\addlinespace " + col if [_n]>=2
	order col
	foreach var of varlist B* Rsq {
		qui tostring `var', force replace format("%5.3f")
	}	
	tempfile tmp
	save "`tmp'"
	
	* Put t-stats below coefficients		
	keep col SE*
	forval i = 1/`num_vars' {
		ren SE`i' B`i'
		}
	append using "`tmp'"
	gsort col -SE1
	replace col = "" if mi(SE1)

	* Rename columns
	local j = 1
	foreach var in `state_vars' {
		if "`var'"=="scale" {
			replace col = "\addlinespace $  iva_{t} $" if [_n]==`j'	 
		}
		else {
			replace col = "\addlinespace $  `var'_{t} $" if [_n]==`j'									
		}
		local j = `j' + 2
	} 
	drop SE*

	* Export to Latex 
	local columns = ""
	local c = "c"
	local total_vars = "`state_vars'" 
	
													
	foreach var in `state_vars' {
		if "`var'"=="scale" {
			local columns = "`columns' & $  iva_{t-1} $"	 
		}
		else {
			local columns = "`columns' & $ `var'_{t-1} $"									
		}
		local c = "`c' c"
	}				
	local columns = "`columns' & $ R^{2} $" 

	* Transition matrix
	listtab * using "Results/Tables/02_06_`vweight'_`end'.tex", ///
		rstyle(tabular) replace ///
		head("\begin{tabular}{l`c'}" ///
		"\midrule" ///
		"`columns' \\" ///
		"\midrule") ///
		foot("\midrule" "\end{tabular}") 
	
		 
****************************************************************************
* 02.07  Compute variance-covariance matrix of residuals  (Table 02_07_`vweight')
**************************************************************************** 		 
	local j = 0
	local varlist `state_vars'  
	foreach var1 in `varlist' {
		local j = `j' + 1
		local k = 0
		foreach var2 in `varlist' {
			local k = `k' + 1
				
			* Weighted average of u_`var1' and u_`var2'   
			u "`resid'", clear
			g one = 1
			tempfile tmp
			save "`tmp'"
			collapse mean1 = u_`var1' mean2 = u_`var2', by(year one)
			collapse mean1 mean2, by(one)
			qui merge 1:m one using "`tmp'", nogen
				
			* Compute weighted variance (covariance) and save locally 
			g one_weight = weight
			qui g cov_`j'_`k' = (u_`var1'-mean1) * (u_`var2'-mean2) * one_weight
			collapse (sum) cov_`j'_`k' one_weight
			local cov_`j'_`k' = cov_`j'_`k'[1] / one_weight[1]

		}
	}  
	* Export to latex
	clear
	svmat B
	local numrows = _N 
	forval i = 1/`numrows' {
		forval j = 1/`numrows' {  
			replace B`i' = `cov_`i'_`j'' if [_n]==`j'
		}
	} 
	local numrows_ = `numrows' + 1
	
	mkmat B*, mat(Sigma) 
	g col = ""
	local i = 1
	foreach var in `state_vars' {
		if "`var'"=="scale" {
			replace col = " $ iva_{t} $" if [_n]==`i'
		}
		else {
			replace col = " $ `var'_{t} $" if [_n]==`i'
		}
		local i = `i' + 1
	}
	 
	* Keep only state variables
	replace col = "\addlinespace " + col if [_n]>=2
	order col
	foreach var of varlist B* {
		qui tostring `var', force replace format("%5.3f")
	} 
	local columns = ""
	local c = "c" 
	foreach var in `state_vars' {
		if "`var'"=="scale" {
			local columns = "`columns' & $ iva_{t} $"
		}
		else {
			local columns = "`columns' & $ `var'_{t} $"
		}
		local c = "`c' c"
	}	
	
	listtab * using "Results/Tables/02_07_`vweight'_`end'.tex", ///
		rstyle(tabular) replace ///
		head("\begin{tabular}{l`c'}" ///
		"\midrule" ///
		"`columns' \\" ///
		"\midrule") ///
		foot("\midrule" "\end{tabular}")

****************************************************************************
* 02.08  Decomposition of return news  (Tables 02_08_`vweight'_`count')
**************************************************************************** 	 

	* Define elementary vectors  
	forvalues i = 1/3 {
		matrix e`i' = J(`num_vars',1,0) 
		matrix e`i'[`i',1] = 1 
	} 
 
	* Generate lambdas 
	matrix lambda =	(e1'*`rho'*B*inv(I(`num_vars')-`rho'*B))'
	matrix lambda2 = (e2'*inv(I(`num_vars')-`rho'*B))'						// For computing roe news explicitly
	matrix lambda3 = (e3'*inv(I(`num_vars')-`rho'*B))'						// For computing scale news explicitly
	matrix lambda2r = (e1 + lambda - lambda3)								// For computing roe news as a residual term
	matrix lambda3r = (e1 + lambda - lambda2)								// For computing scale news as a residual term
	matrix lambda4r = (e1 + lambda - lambda2 -lambda3)								// For computing scale news as a residual term	

	* Generate news variance and covariance
	forvalues count = 1/3 {
		* `count'==1: Back out everything explicitly plus pva as residual
	
		* Variances 
		matrix var_N_DR = lambda'*Sigma*lambda
		matrix var_N_roe = lambda2'*Sigma*lambda2
		matrix var_N_scale = lambda3'*Sigma*lambda3
		matrix var_N_pva = lambda4r'*Sigma*lambda4r
	
		* Covariances	
		matrix cov_N_DR_roe = lambda'*Sigma*lambda2
		matrix cov_N_DR_scale = lambda'*Sigma*lambda3
		matrix cov_N_roe_scale = lambda2'*Sigma*lambda3	
		matrix cov_N_DR_pva = lambda'*Sigma*lambda4r
		matrix cov_N_roe_pva = lambda2'*Sigma*lambda4r
		matrix cov_N_scale_pva = lambda3'*Sigma*lambda4r	
/*	
		* Backing out one of the news as the residual
		if `count'==2 {
			// Back out roe news as the residual
			matrix var_N_roe = lambda2r'*Sigma*lambda2r
			matrix cov_N_DR_roe = lambda'*Sigma*lambda2r
			matrix cov_N_roe_scale = lambda2r'*Sigma*lambda3	
		}
		if `count'==3 {
			// Back out scale news as the residual
			matrix var_N_scale = lambda3r'*Sigma*lambda3r
			matrix cov_N_DR_scale = lambda'*Sigma*lambda3r
			matrix cov_N_roe_scale = lambda2'*Sigma*lambda3r	
		}
*/
		* Variance-covariance matrix	
		local var var_N_DR var_N_roe var_N_scale var_N_pva cov_N_DR_roe cov_N_DR_scale cov_N_roe_scale cov_N_DR_pva cov_N_roe_pva cov_N_scale_pva
		foreach mat in `var' {
			local `mat' = `mat'[1,1]
		}
		matrix totalvar = var_N_DR + var_N_roe + var_N_scale - 2*(cov_N_DR_roe + cov_N_DR_scale - cov_N_roe_scale)
		matrix totalvar_pva = var_N_DR + var_N_roe + var_N_scale + var_N_pva - 2*(cov_N_DR_roe + cov_N_DR_scale - cov_N_roe_scale + cov_N_DR_pva -cov_N_roe_pva - cov_N_scale_pva)
		local totalvar = totalvar[1,1]
		local totalvar_pva = totalvar_pva[1,1]		
		local svarDR = `var_N_DR' / `totalvar'
		local svarroe = `var_N_roe' / `totalvar'
		local svarscale = `var_N_scale' / `totalvar'
		local svarDRroe = -2*`cov_N_DR_roe' / `totalvar'
		local svarDRscale = -2*`cov_N_DR_scale' / `totalvar'
		local svarroescale = 2*`cov_N_roe_scale' / `totalvar'
		
		local svarDR_pva = `var_N_DR' / `totalvar_pva'
		local svarroe_pva = `var_N_roe' / `totalvar_pva'
		local svarscale_pva = `var_N_scale' / `totalvar_pva'
		local svarpva_pva = `var_N_pva' / `totalvar_pva'		
		local svarDRroe_pva = -2*`cov_N_DR_roe' / `totalvar_pva'
		local svarDRscale_pva = -2*`cov_N_DR_scale' / `totalvar_pva'
		local svarroescale_pva = 2*`cov_N_roe_scale' / `totalvar_pva'
		local svarDRpva_pva = -2*`cov_N_DR_pva' / `totalvar_pva'
		local svarroepva_pva = 2*`cov_N_roe_pva' / `totalvar_pva'
		local svarscalepva_pva = 2*`cov_N_scale_pva' / `totalvar_pva'		
		
		local var var_N_DR var_N_roe var_N_scale var_N_pva
		foreach mat in `var' {
			local `mat' = sqrt(`mat'[1,1])
		}
		local cov_N_DR_roe = `cov_N_DR_roe'/(`var_N_DR'*`var_N_roe')
		local cov_N_DR_scale = `cov_N_DR_scale'/(`var_N_DR'*`var_N_scale')
		local cov_N_roe_scale = `cov_N_roe_scale'/(`var_N_roe'*`var_N_scale')
		local cov_N_DR_pva = `cov_N_DR_pva'/(`var_N_DR'*`var_N_pva')
		local cov_N_roe_pva = `cov_N_roe_pva'/(`var_N_roe'*`var_N_pva')
		local cov_N_scale_pva = `cov_N_scale_pva'/(`var_N_scale'*`var_N_pva')		
	
		mat VCM = (`var_N_DR', 0, 0, 0, `svarDR', 0, 0 \ ///
				   `cov_N_DR_roe', `var_N_roe', 0,0, `svarDRroe', `svarroe', 0 \ ///
				   `cov_N_DR_scale', `cov_N_roe_scale', `var_N_scale', 0, `svarDRscale', `svarroescale', `svarscale')
		mat VCM_pva = (`var_N_DR', 0, 0, 0, 0, `svarDR_pva', 0, 0, 0 \ ///
				   `cov_N_DR_roe', `var_N_roe', 0,0,0, `svarDRroe_pva', `svarroe_pva', 0,0 \ ///
				   `cov_N_DR_scale', `cov_N_roe_scale', `var_N_scale', 0,0, `svarDRscale_pva', `svarroescale_pva', `svarscale_pva',0 \ ///		   
 				   `cov_N_DR_pva', `cov_N_roe_pva', `cov_N_scale_pva',`var_N_pva' , 0, `svarDRpva_pva', `svarroepva_pva', `svarscalepva_pva', `svarpva_pva')
				   
		* Export to Latex 
		clear
		svmat VCM
		g col = ""
		g row = 0
		local i = 1
		foreach var in "N_{DR}" "N_{roe}" "N_{iva}" {
				replace col = "$ `var' $" if [_n]==`i'
				replace row = `i'-0.5 if [_n]==`i'
				local i = `i' + 1
				}
		replace col = "\addlinespace " + col if [_n]>=2
		order col 

*			append using "Output/news_decomp_stderr_C.dta"	

		forval i= 1/7 {
			qui tostring VCM`i', force replace format("%5.3f")
*				qui replace VCM`i' = "(" + VCM`i' + ")" if [_n]>3
			qui replace VCM`i' = "" if VCM`i' == "0.000" 
*				qui replace VCM`i' = "" if VCM`i' == "(0.000)"
		}

		sort row
		drop row

		* Output to latex
		listtab * using "Results/Tables/02_08_`vweight'_`count'_`end'.tex", ///
			rstyle(tabular) replace ///
			head("\begin{tabular}{lccccccc}" ///
			"\midrule" ///
			"& \multicolumn{3}{c}{$\sigma$ (diag), $\rho$ (off-diag)} & & \multicolumn{3}{c}{Contribution to $\sigma^2_{N_r}$} \\" ///
			"\cmidrule(lr){2-4} \cmidrule(lr){6-8}" /// 
				" & $ N_{DR} $ & $ N_{roe} $ & $ N_{iva} $ & & $ -N_{DR} $ & $ N_{roe} $ & $ N_{iva} $\\" ///
			"\midrule") ///
			foot("\midrule" "\end{tabular}")
	} 	

		* Export to Latex with pva 
		clear
		svmat VCM_pva
		g col = ""
		g row = 0
		local i = 1
		foreach var in "N_{DR}" "N_{roe}" "N_{iva}" "N_{pva}" {
				replace col = "$ `var' $" if [_n]==`i'
				replace row = `i'-0.5 if [_n]==`i'
				local i = `i' + 1
				}
		replace col = "\addlinespace " + col if [_n]>=2
		order col 

*			append using "Output/news_decomp_stderr_C.dta"	

		forval i= 1/9 {
			qui tostring VCM_pva`i', force replace format("%5.3f")
*				qui replace VCM`i' = "(" + VCM`i' + ")" if [_n]>3
			qui replace VCM_pva`i' = "" if VCM_pva`i' == "0.000" 
*				qui replace VCM`i' = "" if VCM`i' == "(0.000)"
		}

		sort row
		drop row

		* Output to latex
		listtab * using "Results/Tables/02_08_`vweight'_pva_`end'.tex", ///
			rstyle(tabular) replace ///
			head("\begin{tabular}{lccccccccc}" ///
			"\midrule" ///
			"& \multicolumn{4}{c}{$\sigma$ (diag), $\rho$ (off-diag)} & & \multicolumn{4}{c}{Contribution to $\sigma^2_{N_r}$} \\" ///
			"\cmidrule(lr){2-4} \cmidrule(lr){6-8}" /// 
				" & $ N_{DR} $ & $ N_{roe} $ & $ N_{iva} $ & $ N_{pva} $ & & $ -N_{DR} $ & $ N_{roe} $ & $ N_{iva} $ & $ N_{pva} $ \\" ///
			"\midrule") ///
			foot("\midrule" "\end{tabular}")
	
****************************************************************************
* 02.09  R2 analysis   (Table 02_09_`vweight')
****************************************************************************  

	use "`resid'", clear
	g N_return = u_r
	g N_DR=0
	g N_roe=0
	g N_scale=0
	local k=1
	foreach var in `state_vars'{
		replace N_DR = N_DR + lambda[`k',1]*u_`var'
		replace N_roe = N_roe + lambda2[`k',1]*u_`var'
		replace N_scale = N_scale + lambda3[`k',1]*u_`var'
		local k=`k'+1
	}
	g N_CFV = N_return + N_DR // Cash-flow news	
	g N_roe_impl = N_CFV - N_scale // error from N_R on N_DR and N_scale with fixed coefficients
	g N_scale_impl = N_CFV - N_roe // error from N_R on N_DR and N_roe with fixed coefficients
	g N_DR_impl = N_return - N_roe - N_scale // error from N_R on N_roe and N_scale with fixed coefficients
	g N_error_impl = N_CFV - N_roe - N_scale // error from N_R on N_DR, N_roe and N_scale with fixed coefficients
	
	* Label News
	label var 	N_DR			" $ N_{DR} $ "
	label var	N_roe			" $ N_{roe} $ "
	label var 	N_scale			" $ N_{iva} $ "
	
	* Save temporarily
	tempfile tmp
	save "`tmp'"

	/*
	* News term R2  
	eststo, title("$ N_{r} $"): reghdfe N_return N_DR N_roe [aw=weight], nocon absorb(year) cl(year firm)
	eststo, title("$ N_{r} $"): reghdfe N_return N_DR N_roe N_scale [aw=weight], nocon absorb(year) cl(year firm)
	eststo, title("$ N_{CF} $"): reghdfe N_CFV N_roe [aw=weight], nocon absorb(year) cl(year firm)
	eststo, title("$ N_{CF} $"): reghdfe N_CFV N_roe N_scale [aw=weight], nocon absorb(year) cl(year firm)
	 
	* Generate a table
	esttab using "Output/Tables/02_09_`vweight'.tex", replace fragment nostar r2 b(3) se /*t(1)*/ ///
		width(\textwidth) mlabels(,titles)  nonumbers nonotes nolines label substitute(\_ _) ///
		prehead("\begin{tabular}{lcccc} \\ \midrule") posthead(\hline) prefoot(\hline) 
	eststo clear

	* Compute R2 with coefficients fixed at -1 / 1	 
	collapse (sd) N_return N_CF = N_CFV N_scale_impl N_roe_impl N_DR_impl N_error_impl [aw=weight]
	foreach var in N_return N_CF N_scale_impl N_roe_impl N_DR_impl N_error_impl {
		replace `var' = `var'^2
	}
	g R2_scale_r = 1 - N_scale_impl/N_return
*		g R2_roe_r = 1 - N_roe_impl/N_return
*		g R2_DR_r = 1 - N_DR_impl/N_return
	g R2_error_r = 1 - N_error_impl/N_return
	g R2_scale_CF = 1 - N_scale_impl/N_CF
	g R2_error_CF = 1 - N_error_impl/N_CF
*		g R2_roe_CF = 1 - N_roe_impl/N_CF

	local i = 1
	foreach var in R2_scale_r /*R2_roe_r R2_DR_r*/ R2_error_r R2_scale_CF R2_error_CF {
		rename `var' r2_`i'
		local i = `i' + 1
	}
	keep r2_*
	g row = 1.5
*		append using "Output/rsq.dta"	  
	sort row
	drop row 
	foreach var of varlist _all {
		qui tostring `var', force replace format("%5.3f")
*			qui replace `var' = "(" + `var' + ")" if [_n]!=2
		} 
	g col = "" 
	replace col = "$ \tilde{R}^{2}$" /*if [_n]==2*/
	order col	
	* Append to existing table 02_09_`vweight'
	listtab *, appendto("Output/Tables/02_09_`vweight'.tex") rstyle(tabular) ///
		foot("\midrule" "\end{tabular}") 	
	*/	
****************************************************************************
* 02.10  R2 contribution by characteristic  (Table 02_10_`vweight'_`count') (1 is high, 2 is low characteristic table)
**************************************************************************** 				 
/*
	local count = 1
	foreach sign in ">=9" "<=2" {
		use "`tmp'", clear
		local c = ""
		foreach char in `sort_chars' {
			eststo, title("$ N_{CF} $"): reghdfe N_CFV N_roe [aw=weight] if bin_`char' `sign', nocon absorb(year) cl(year firm)
			eststo, title("$ N_{CF} $"): reghdfe N_CFV N_roe N_scale [aw=weight] if bin_`char' `sign', nocon absorb(year) cl(year firm)
			local c = "`c' c c"
		}		
		* Generate a table
		esttab using "Output/Tables/02_10_`vweight'_`count'.tex", replace fragment nostar r2 b(3) se /*t(1)*/ ///
			width(\textwidth) mlabels(,titles)  nonumbers nonotes nolines label substitute(\_ _) ///
			prehead("\begin{tabular}{l`c'} \\ \midrule") posthead(\hline) prefoot(\hline) 
		eststo clear
		 
		* Compute R2 with coefficients fixed at -1 / 1	
		local Rsq = "$ \tilde{R}^{2}$"
		foreach char in `sort_chars' {
			use "`tmp'", clear
			collapse (sd) N_CF = N_CFV N_scale_impl N_error_impl [aw=weight] if bin_`char' `sign'
			foreach var in N_CF N_scale_impl N_error_impl {
				replace `var' = `var'^2
			} 
			g R2_error_CF = 1 - N_error_impl/N_CF
			g R2_scale_CF = 1 - N_scale_impl/N_CF 
		
			* Save R2 as local variables  
			local rsq = R2_scale_CF[1]
			local rsq: di %04.3f `rsq'
			local Rsq = "`Rsq' & `rsq'" 
			local rsq = R2_error_CF[1]
			local rsq: di %04.3f `rsq'
			local Rsq = "`Rsq' & `rsq'" 
			display "`Rsq'" 
		}
		clear
		gen R2theory = ""
		set obs `=_N+1'
		replace R2theory = "`Rsq'"
		listtab R2theory, appendto("Output/Tables/02_10_`vweight'_`count'.tex") rstyle(tabular) ///
			foot("\midrule" "\end{tabular}")
		// Read the exported table and change constrained indicator
		infix str l 1-500 using ///
			"Output/Tables/02_10_`vweight'_`count'.tex", clear
			replace l = subinstr(l,"&  $ N_{CF} $&  $ N_{CF} $&  $ N_{CF} $&  $ N_{CF} $&  $ N_{CF} $&  $ N_{CF} $&  $ N_{CF} $&  $ N_{CF} $","`subsample_header`count''",.)  
			// Drop number of observations & standard errors on news terms
			drop if inlist([_n],5,8,10)
			if `count'==2 {
				drop if inlist([_n],11)
			}
			// Replace
			outfile using ///
				"Output/Tables/02_10_`vweight'_`count'.tex", ///
				noquote replace 
			
			if `count'==2 { 
/*				u "Data/created/bootstrap R2 CI.dta", clear
				gen R2_gap1_loA = R2_gap1_lo
				gen R2_gap4_loA = R2_gap4_lo				
				replace R2_gap1_lo = -R2_gap1_hi
				replace R2_gap4_lo = -R2_gap4_hi
				replace R2_gap1_hi = -R2_gap1_loA
				replace R2_gap4_hi = -R2_gap4_loA	
				forval i=1/4 {
					qui tostring R2_gap`i'_lo, force replace format("%5.3f")
					qui tostring R2_gap`i'_hi, force replace format("%5.3f")
					gen CI`i'= "\multicolumn{2}{c}{[" + R2_gap`i'_lo + ", " + R2_gap`i'_hi + "]}"
				}
				drop R2*
				gen CI = "CI($\Delta R^2$)" if row==1
				replace CI = "CI($\Delta\tilde{R}^2$)" if mi(CI)
*/
				u "Data/created/bootstrap R2 p_recentered.dta", clear
				forval i=1/5 {
					qui tostring P`i', force replace format("%5.3f")
					gen CI`i'= "[" + P`i' + "] & "
				}
				gen CI = "$p(\Delta R^2)$" if row==1
				replace CI = "$p(\Delta\tilde{R}^2)$" if mi(CI)
				drop row
				order CI
				listtab CI*, delimiter(&) appendto("Output/Tables/02_10_`vweight'_`count'.tex") rstyle(tabular) ///
				foot("\midrule" "\end{tabular}")

				}	
			* Update count
			local count = `count' + 1
	}
	local vweight = 0
	// Also generate a version of table 02_10_`vweight'_2 without p values
	infix str l 1-500 using "Output/Tables/02_10_`vweight'_2.tex", clear	
		drop if inlist([_n],10,11,12) 
	// Replace  
	outfile using "Output/Tables/02_10_`vweight'_3.tex", noquote replace
*/
****************************************************************************
* 02.10.2  Variance decomposition by characteristic, with full-sample transition matrix (1 is high, 2 is low characteristic table)
**************************************************************************** 				 
if `subsamples'==1 {
	local count = 1
	foreach sign in ">=9" "<=2" {
		foreach char in `sort_chars' {
			use "`tmp'", clear
			keep if bin_`char' `sign'
			foreach news1 in "DR" "roe" "scale"{
				foreach news2 in "DR" "roe" "scale"{
					correlate N_`news1' N_`news2', cov
					matrix cov_N_`news1'_`news2' = r(cov_12)
					}
				matrix var_N_`news1' = cov_N_`news1'_`news1'
				}
			* Variance-covariance matrix	
			local var var_N_DR var_N_roe var_N_scale cov_N_DR_roe cov_N_DR_scale cov_N_roe_scale
			foreach mat in `var' {
				local `mat' = `mat'[1,1]
			}
			matrix totalvar = var_N_DR + var_N_roe + var_N_scale - 2*(cov_N_DR_roe + cov_N_DR_scale - cov_N_roe_scale)
			local totalvar = totalvar[1,1]
			local svarDR = `var_N_DR' / `totalvar'
			local svarroe = `var_N_roe' / `totalvar'
			local svarscale = `var_N_scale' / `totalvar'
			local svarDRroe = -2*`cov_N_DR_roe' / `totalvar'
			local svarDRscale = -2*`cov_N_DR_scale' / `totalvar'
			local svarroescale = 2*`cov_N_roe_scale' / `totalvar'			
			local var var_N_DR var_N_roe var_N_scale
			foreach mat in `var' {
				local `mat' = sqrt(`mat'[1,1])
			}
			local cov_N_DR_roe = `cov_N_DR_roe'/(`var_N_DR'*`var_N_roe')
			local cov_N_DR_scale = `cov_N_DR_scale'/(`var_N_DR'*`var_N_scale')
			local cov_N_roe_scale = `cov_N_roe_scale'/(`var_N_roe'*`var_N_scale')	
		
			mat VCM_sub = (`var_N_DR', 0, 0, 0, `svarDR', 0, 0 \ ///
					   `cov_N_DR_roe', `var_N_roe', 0,0, `svarDRroe', `svarroe', 0 \ ///
					   `cov_N_DR_scale', `cov_N_roe_scale', `var_N_scale', 0, `svarDRscale', `svarroescale', `svarscale')
	 
			* Export to Latex 
			clear
			svmat VCM_sub
			g col = ""
			g row = 0
			local i = 1
			foreach var in "N_{DR}" "N_{roe}" "N_{iva}" {
					replace col = "$ `var' $" if [_n]==`i'
					replace row = `i'-0.5 if [_n]==`i'
					local i = `i' + 1
					}
			replace col = "\addlinespace " + col if [_n]>=2
			order col 

			forval i= 1/7 {
				qui tostring VCM_sub`i', force replace format("%5.3f")
				qui replace VCM_sub`i' = "" if VCM_sub`i' == "0.000" 
			}

			sort row
			drop row

			* Output to latex
			listtab * using "Results/Tables/02_10_`char'_`count'_`end'.tex", ///
				rstyle(tabular) replace ///
				head("\begin{tabular}{lccccccc}" ///
				"\midrule" ///
				"& \multicolumn{3}{c}{$\sigma$ (diag), $\rho$ (off-diag)} & & \multicolumn{3}{c}{Contribution to $\sigma^2_{N_r}$} \\" ///
				"\cmidrule(lr){2-4} \cmidrule(lr){6-8}" /// 
					" & $ N_{DR} $ & $ N_{roe} $ & $ N_{iva} $ & & $ -N_{DR} $ & $ N_{roe} $ & $ N_{iva} $\\" ///
				"\midrule") ///
				foot("\midrule" "\end{tabular}")
		}
		* Update count
		local count = `count' + 1
	}
}

****************************************************************************
* 02.11  Decompose the cross-sectional var(mb)   (Table 02_11_`vweight')
**************************************************************************** 	 
	use "`resid'", clear	 
	if `vweight2'==1{
		replace weight = valueweight
	}
	matrix expected = inv(I(`num_vars')-`rho'*B) * B
		mat list expected
		
	matrix expected_r = expected[1, 1..`num_vars']
	matrix expected_roe = expected[2, 1..`num_vars'] 
	matrix expected_scale = expected[3, 1..`num_vars']
	foreach var in r roe scale {
		svmat expected_`var'
		forvalues i = 1/`num_vars' {
			replace expected_`var'`i' = expected_`var'`i'[_n-1] if [_n]>1
		}
		gen exp_`var' = r * expected_`var'1
		local i = 2
		foreach var2 in roe scale bm {
			replace exp_`var' = exp_`var' + `var2' * expected_`var'`i'
			local i = `i' + 1
		}
		drop expected_`var'*
	} 
	
	* Output expected_r, expected_roe, expected_scale for use in other research
	tempfile tmp
	save "`tmp'"	
	keep year gvkey exp_r-exp_scale
	save "Results/Other/expected_cf", replace
	use "`tmp'", clear
	
	* Cross-sectional decomposition
	local i = 1
	replace exp_r = -1 * exp_r
	foreach var in r roe scale {
		label var mb "$ mb_t $"
		reghdfe exp_`var' mb [aweight = weight], nocon absorb(year) cl(year firm) // if year<=`end'-15  
			* Save results	  
			eststo C`i'
			local i = `i' + 1	
	}
	// Tentative table in latex
	esttab C1 C2 C3 using ///
		"Results/Tables/02_11_`vweight'_`vweight2'_`end'.tex", ///
		 b(3) se /*t(1)*/  margin /*booktabs*/ nomtitles ///
		 nostar label ///
		 mgroups("$ -\sum_{j=1}^{J}\rho^{j-1} E_t\left[r_{t+j}\right] $" ///
			"$ \sum_{j=1}^{J}\rho^{j-1} E_t\left[roe_{t+j}\right] $" ///
			"$ \sum_{j=1}^{J}\rho^{j-1} E_t\left[iva_{t+j}\right] $", pattern(1 1 1)                   ///
		 prefix(\multicolumn{@span}{c}{) suffix(})   ///
		 span erepeat(\cmidrule(lr){@span}))    /// 
		mlabels(none) collabels(none) stats(N, fmt(%5.0fc) ///
		labels("Observations")) /// 
		nonotes replace constant substitute(\_ _) 	
		
	preserve	
		
	infix str l 1-500 using "Results/Tables/02_11_`vweight'_`vweight2'_`end'.tex", clear	 
		replace l = "\hline" if l=="\hline\hline"
		replace l = " $ J$ &\multicolumn{1}{c}{$ -\sum_{j=1}^{J}\rho^{j-1} E_t\left[r_{t+j}\right] $}&\multicolumn{1}{c}{$ \sum_{j=1}^{J}\rho^{j-1} E_t\left[roe_{t+j}\right] $}&\multicolumn{1}{c}{$ \sum_{j=1}^{J}\rho^{j-1} E_t\left[iva_{t+j}\right] $}\\" if [_n]==3
		drop if inlist([_n],4,8,9)
		replace l = subinstr(l, "mb_t", "\infty", .) if [_n]==5  
	// Replace 
	outfile using "Results/Tables/02_11_`vweight'_`vweight2'_`end'.tex", noquote replace
	 
	eststo clear

	restore
	
	// Generate coefficients for Table C6 
	reghdfe exp_r mb exp_roe exp_scale [aweight = weight], nocon absorb(year) cl(year firm) 				
		matrix B = e(b)
		matrix V = vecdiag(e(V))
		matmap V V, map(sqrt(@))
		matrix Rsq = e(r2)
		matrix N = e(N)
		clear 
		svmat B
		svmat V
		svmat Rsq
		svmat N
	export excel "Results/Other/02_11_coeff_stderr.xlsx", replace
		

****************************************************************************
* 02.12  MVE portfolio analysis  (Tables 02_12_`vweight'_1 through 3)
**************************************************************************** 
	* Compute median and the number of firms in each year
	use "`resid'", clear
	collapse (median) `all_vars', by(year)
		
	* Rename median values
	foreach var in `all_vars' {	
		ren `var' median`var'														
	}	
	merge 1:m year using "`resid'", nogen
	drop u_*
	
	* Define regressors and set of variables including medians	 
	local regressors = ""
	local regressors2 = ""															// Median regressors	
	local median_vars = ""   														// Set of medians
	foreach var in `state_vars' {
		forval i = 1/`lags' { 
			local regressors = "`regressors' `var'_L`i'"
			local regressors2 = "`regressors2' median`var'_L`i'"					// Regrssors with restriction
			local median_vars = "`median_vars' median`var'"
		}
	}	 
	local regressors = "`regressors' `regressors2'"
	
	* Count the number of median values
	local num_vars = 0
	local medians = 0
	foreach var in `state_vars' {
		local num_vars = `num_vars' + 2									// +2 since we are adding median values for each state variable
		local medians = `medians' + 1
	}

	* Define matrix to store the results	
	matrix B = J(1,`num_vars',.)													// Coefficient matrix
	matrix SE = J(1,`num_vars',.)													// Standard error matrix
	matrix Rsq = J(1,1,.)

	* Regress state variables on lagged values and lagged medians
	foreach var of varlist `state_vars' {
		
		* Regression 
		qui reghdfe `var' `regressors' [pweight = weight], nocon noabsorb cl(year firm) resid
					
		* Compute error terms	 
		predict u_`var', residual
					
		* Store coefficients and standard errors 
		matrix B = (B \ (e(b)))
		matrix V = vecdiag(e(V))
		matmap V V, map(sqrt(@))
		matrix SE = (SE \ V)	
		matrix Rsq = (Rsq \ (e(r2)))						
	}
	
	* Regress medians with restrictions (an individual firm does not affect medians)
	 foreach var of varlist `median_vars' {
		
		* Regression 
		qui reghdfe `var' `regressors2', nocon noabsorb cl(year) resid
					
		* Compute error terms	 
		predict u_`var', residual

		* Store coefficients and standard errors 
		matrix B = [B \ [J(1,`medians',0),(e(b))]]
		matrix V = [J(1,`medians',0),vecdiag(e(V))]
		matmap V V, map(sqrt(@))
		matrix SE = (SE \ V)	
		matrix Rsq = (Rsq \ (e(r2)))				
	}
	matrix B = (B, Rsq)			

	* Decomposition of return news (back out all news explicitly)
	matrix e1 = [1\0\0\0\J(`num_vars'-4,1,0)]
	matrix e2 = [0\1\0\0\J(`num_vars'-4,1,0)]
	matrix e3 = [0\0\1\0\J(`num_vars'-4,1,0)] 
	
	* Generate lambdas 
	local num_vars2 = `num_vars' + 1
	matrix B = B[2..`num_vars2',1..`num_vars'] 
	matrix lambda =	(e1'*`rho'*B*inv(I(`num_vars')-`rho'*B))'
	matrix lambda2 = (e2'*inv(I(`num_vars')-`rho'*B))'					// For computing roe news explicitly
	matrix lambda3 = (e3'*inv(I(`num_vars')-`rho'*B))'					// For computing scale news explicitly

	* Compute news term
	g N_return = u_r
	g N_DR=0
	g N_scale=0
	g N_roe=0 
	local k=1
	foreach var in `state_vars' `median_vars'{
		replace N_DR = N_DR + lambda[`k',1]*u_`var'
		replace N_roe = N_roe + lambda2[`k',1]*u_`var'
		replace N_scale = N_scale + lambda3[`k',1]*u_`var'
		local k=`k'+1
	} 
	g N_CFV = N_return + N_DR
	
	* Save data
	tempfile NewsIntermediate
	save "`NewsIntermediate'"
	 	
	* Compute News terms of anomaly portfolios
	local anomalies = "bm_L1 ME_L1 ROE_L1 assetgrowth_L1 mom_L1 "
	local charloop = 1
	foreach char in `anomalies' {
		u "`NewsIntermediate'", clear
		* Subsample	high
		keep if bin_`char'> 8
		collapse (mean) N_returnH = N_return N_DRH = N_DR N_roeH = N_roe N_scaleH = N_scale [pweight=ME_L1], by(year)
		tempfile NewsHigh
		save "`NewsHigh'", replace

		* Subsample	low
		u "`NewsIntermediate'", clear
		keep if bin_`char'< 3
		collapse (mean) N_returnL = N_return N_DRL = N_DR N_roeL = N_roe N_scaleL = N_scale [pweight=ME_L1], by(year)
		merge 1:1 year using "`NewsHigh'", nogen
		
		foreach var in N_return N_DR N_roe N_scale {
			g `var'`char' = `var'H - `var'L
			drop `var'H `var'L
			}
		tempfile NewsAnomalies`char'
		save "`NewsAnomalies`char''" 
		local N_vars = "N_DR`char' N_roe`char' N_scale`char'"
		qui corr `N_vars', cov
		mat A = r(C)
		qui corr `N_vars'
		mat B = r(C)
		local totalvar = A[1,1] + A[2,2] + A[3,3] - 2*(A[2,1] + A[3,1] - A[3,2])
		mat VCM = (A[1,1]^(0.5), 0, 0, 0 , A[1,1]/`totalvar', 0, 0 \ ///
				   B[2,1], A[2,2]^(0.5), 0, 0, -2*A[2,1]/`totalvar', A[2,2]/`totalvar', 0 \ ///
				   B[3,1], B[3,2], A[3,3]^(0.5), 0, -2*A[3,1]/`totalvar', 2*A[3,2]/`totalvar', A[3,3]/`totalvar') 
		clear
		svmat VCM
		g col = ""
		local i = 1
		foreach var in "N_{DR}" "N_{roe}" "N_{iva}" {
				replace col = "$ `var' $" if [_n]==`i'
				local i = `i' + 1
				}
		replace col = "\addlinespace " + col if [_n]>=2
		order col
		foreach var of varlist VCM* {
			qui tostring `var', force replace format("%5.3f")
			} 
		forval i= 1/7{
			replace VCM`i' = "" if VCM`i' == "0.000" 
		}
		if `charloop'==1 {
		listtab * using "Results/Tables/02_12_`vweight'_1_`char'_`end'.tex", /// ///
			rstyle(tabular) replace ///
			head("\begin{tabular}{lccccccc}" ///
			"\midrule" ///
			"& \multicolumn{3}{c}{$\sigma$ (diag), $\rho$ (off-diag)} & & \multicolumn{3}{c}{Contribution to $\sigma^2_{N_r}$} \\" ///
		"\cmidrule(lr){2-4} \cmidrule(lr){6-8}" /// 
		" & $ N_{DR} $ & $ N_{roe} $ & $ N_{iva} $ & & $ -N_{DR} $ & $ N_{roe} $ & $ N_{iva} $ \\" ///
			"\midrule") ///
			foot("\midrule" "\end{tabular}")
		}
		else {
		listtab * using "Results/Tables/02_12_`vweight'_1_`char'_`end'.tex", /// ///
			rstyle(tabular) replace ///
			head("\begin{tabular}{lccccccc}" ///
			"\midrule") ///
			foot("\midrule" "\end{tabular}")			
		}
		local charloop = `charloop' + 1		
		u "`NewsAnomalies`char''", clear	
	}
	tempfile NewsAnomalies
	save "`NewsAnomalies'"

	foreach char in bm_L1 ME_L1 ROE_L1 assetgrowth_L1 {
		merge 1:1 year using "`NewsAnomalies`char''", nogen
		save "`NewsAnomalies'", replace
	}
	u "`NewsIntermediate'", clear
	collapse (mean) N_returnM = N_return N_DRM = N_DR N_roeM = N_roe N_scaleM = N_scale [pweight=ME_L1], by(year)
	merge 1:1 year using "`NewsAnomalies'", nogen
	local N_vars = "N_DR N_roe N_scale"
	foreach var in `N_vars' {
		g `var'MVE = 0.8*`var'M - 0.21*`var'bm_L1 -0.35*`var'ME_L1 + 0.73*`var'ROE_L1 -1.87*`var'assetgrowth_L1 + 1.35*`var'mom_L1
		g `var'MVEex = 0.06*`var'bm_L1 -0.8*`var'ME_L1 + 0.66*`var'ROE_L1 -1.55*`var'assetgrowth_L1 + 1.44*`var'mom_L1
	}
	save "`NewsIntermediate'", replace	
	local N_varsMVE = "N_DRMVE N_roeMVE N_scaleMVE"

	qui corr `N_varsMVE', cov
	mat A = r(C)
	qui corr `N_varsMVE'
	mat B = r(C)
	local totalvar = A[1,1] + A[2,2] + A[3,3] - 2*(A[2,1] + A[3,1] - A[3,2])
	mat VCM = (A[1,1]^(0.5), 0, 0, 0 , A[1,1]/`totalvar', 0, 0 \ ///
			   B[2,1], A[2,2]^(0.5), 0, 0, -2*A[2,1]/`totalvar', A[2,2]/`totalvar', 0 \ ///
			   B[3,1], B[3,2], A[3,3]^(0.5), 0, -2*A[3,1]/`totalvar', 2*A[3,2]/`totalvar', A[3,3]/`totalvar') 
	clear
	svmat VCM
	g col = ""
	local i = 1
	foreach var in "N_{DR}" "N_{roe}" "N_{iva}" {
			replace col = "$ `var' $" if [_n]==`i'
			local i = `i' + 1
		}
	replace col = "\addlinespace " + col if [_n]>=2
	order col
	foreach var of varlist VCM* {
		qui tostring `var', force replace format("%5.3f")
	} 
	forval i= 1/7{
		replace VCM`i' = "" if VCM`i' == "0.000" 
	} 
	listtab * using "Results/Tables/02_12_`vweight'_2_`end'.tex", /// ///
		rstyle(tabular) replace ///
		head("\begin{tabular}{lccccccc}" ///
		"\midrule") ///
		foot("\midrule" "\end{tabular}")	 
	u "`NewsIntermediate'", clear
	local N_varsMVEex = "N_DRMVEex N_roeMVEex N_scaleMVEex"

	qui corr `N_varsMVEex', cov
	mat D = r(C)
	qui corr `N_varsMVEex'
	mat E = r(C)
	local totalvar = D[1,1] + D[2,2] + D[3,3] - 2*(D[2,1] + D[3,1] - D[3,2])
	mat VCMex = (D[1,1]^(0.5), 0, 0, 0 , D[1,1]/`totalvar', 0, 0 \ ///
			   E[2,1], D[2,2]^(0.5), 0, 0, -2*D[2,1]/`totalvar', D[2,2]/`totalvar', 0 \ ///
			   E[3,1], E[3,2], D[3,3]^(0.5), 0, -2*D[3,1]/`totalvar', 2*D[3,2]/`totalvar', D[3,3]/`totalvar')
	clear		   
	svmat VCMex
	g col = ""
	local i = 1
	foreach var in "N_{DR}" "N_{roe}" "N_{iva}" {
			replace col = "$ `var' $" if [_n]==`i'
			local i = `i' + 1
	}
	replace col = "\addlinespace " + col if [_n]>=2
	order col
	foreach var of varlist VCMex* {
		qui tostring `var', force replace format("%5.3f")
	}		
	forval i= 1/7 {
		replace VCMex`i' = "" if VCMex`i' == "0.000" 
	}
	listtab * using "Results/Tables/02_12_`vweight'_3_`end'.tex", /// ///
		rstyle(tabular) replace ///
		head("\begin{tabular}{lccccccc}" ///
		"\midrule") ///
		foot("\midrule" "\end{tabular}") 
	
	
********************************************************************************
* 02.13  Repeat VAR analysis for spec 2  
******************************************************************************** 
use "`data'", clear
local spec = 2
local lags = 4

* Weight on each observation ((#firms in the year)^(-1) * possible value weight)
tempfile tmp
save "`tmp'" 
collapse (count) firms = mb (sum) ME_L1, by(year)
ren ME_L1 sumME_L1
merge 1:m year using "`tmp'", nogen
g weight = 1/firms
if `vweight'==1 {
	replace weight = ME_L1/sumME_L1
}
gen valueweight = ME_L1/sumME_L1

* Choose state variables	
local state_vars = "`basic_vars' `extra_vars'"

* Choose regressors
local regressors = "" 
forval i = 1/`lags' {
	foreach var in `state_vars' {
		if "`var'"=="bm" & `i'<=1 {
			local regressors = "`regressors' `var'_L`i'" 
		}
		else if inlist("`var'","roe","scale","sgrowth","gamma") & `i'<=2 {
			local regressors = "`regressors' `var'_L`i'" 
		}
		else if inlist("`var'","r","assetgrowth") {
			local regressors = "`regressors' `var'_L`i'"  
		}
	}
}	 

* Count the number of regressors
local num_vars = 0
foreach var in `regressors'{
	local num_vars = `num_vars' + 1
	}

* Choose all variables including lags
local all_vars = "`state_vars'" 
local all_vars = "`all_vars' `regressors'"

* Define matrix to store the results
matrix B = J(1,`num_vars',.)							// Coefficient matrix
matrix SE = J(1,`num_vars',.)							// Standard error matrix
matrix Rsq = J(1,1,.)

****************************************************************************
* 02.14  Regress state variables (WLS)  (Table 02_14_`vweight')
**************************************************************************** 
	
local i = 1
foreach var of varlist `state_vars' {
	 
	* Regression
	reghdfe `var' `regressors' [aweight = weight], nocon absorb(year) cl(year firm) resid
		 
	* Compute error terms				
	predict u_`var', residual
				 
	* Store coefficients and standard errors						
	matrix B = (B \ (e(b)))
	matrix V = vecdiag(e(V))
	matmap V V, map(sqrt(@))
	matrix SE = (SE \ V)
	matrix Rsq = (Rsq \ (e(r2_within)))
 
	local i = `i' + 1	
} 
 
* Remove first row
mat B = B[2...,1...]
mat Rsq = Rsq[2...,1...]
mat SE = SE[2...,1...]
  
* Save data on residuals
tempfile resid
save "`resid'"

* Create a table of coefficients
* Import stored coefficients and standard errors 
clear 
svmat B
svmat Rsq
svmat SE

ren Rsq* Rsq
 
* Standard erros  
forval i = 1/`num_vars'{ 
	qui tostring SE`i', force replace format("%5.3f")
	qui replace SE`i' = "(" + SE`i' + ")"
} 

* Make a column for variable names	
g col = ""
local j = 1
foreach var in `state_vars' {
	if "`var'"=="scale" {
		replace col = "$ 0`j'iva_{t} $" if [_n]==`j'
	}
	else {
		replace col = "$ 0`j'`var'_{t} $" if [_n]==`j'
	}
	local j = `j' + 1
}  
replace col = "\addlinespace " + col if [_n]>=2
order col
foreach var of varlist B* Rsq {
	qui tostring `var', force replace format("%5.3f")
}	
tempfile tmp
save "`tmp'"

* Put t-stats below coefficients		
keep col SE*
forval i = 1/`num_vars' {
	ren SE`i' B`i'
	}
append using "`tmp'"
gsort col -SE1
replace col = "" if mi(SE1)

* Rename columns
local j = 1
foreach var in `state_vars' {
	if "`var'"=="scale" {
		replace col = "\addlinespace $  iva_{t} $" if [_n]==`j'	 
	}
	else if "`var'"=="assetgrowth"{
		replace col = "\addlinespace $  ag_{t} $" if [_n]==`j'	
	}
	else {
		replace col = "\addlinespace $  `var'_{t} $" if [_n]==`j'									
	}
	local j = `j' + 2
} 
drop SE*

* Export to Latex 
local c = "c"								
foreach var in `regressors' {
	local c = "`c' c"
}				

* Transition matrix 
listtab * using "Results/Tables/02_14_`vweight'_`end'.tex", ///
	rstyle(tabular) replace ///
	head("\begin{tabular}{l`c'}" ///
	"\midrule" ///
	"& $ r_{t-1} $ & $ roe_{t-1} $ & $ iva_{t-1} $ & $ bm_{t-1} $ & $ ag_{t-1} $ & $ r_{t-2} $ & $ roe_{t-2} $ & $ iva_{t-2} $ & $ ag_{t-3} $ & $ r_{t-3} $ & $ ag_{t-3} $ & $ r_{t-4} $ & $ ag_{t-4} $ & $ R^2 $ \\" ///
	"\midrule") ///
	foot("\midrule" "\end{tabular}") 
		 
****************************************************************************
* 02.15  Compute variance-covariance matrix of residuals  (Table 02_15_`vweight')
**************************************************************************** 		 
local j = 0
local varlist `state_vars'  
foreach var1 in `varlist' {
	local j = `j' + 1
	local k = 0
	foreach var2 in `varlist' {
		local k = `k' + 1
			
		* Weighted average of u_`var1' and u_`var2'   
		u "`resid'", clear
		g one = 1
		tempfile tmp
		save "`tmp'"
		collapse mean1 = u_`var1' mean2 = u_`var2', by(year one)
		collapse mean1 mean2, by(one)
		qui merge 1:m one using "`tmp'", nogen
			
		* Compute weighted variance (covariance) and save locally 
		g one_weight = weight
		qui g cov_`j'_`k' = (u_`var1'-mean1) * (u_`var2'-mean2) * one_weight
		collapse (sum) cov_`j'_`k' one_weight
		local cov_`j'_`k' = cov_`j'_`k'[1] / one_weight[1]

	}
}  

* Export to latex
clear
svmat B
local numrows = _N 
forval i = 1/`numrows' {
	forval j = 1/`numrows' {  
		replace B`i' = `cov_`i'_`j'' if [_n]==`j'
	}
} 
local numrows_ = `numrows' + 1
drop B`numrows_'-B`num_vars' 

mkmat B*, mat(Sigma) 
g col = ""
local i = 1
foreach var in `state_vars' {
	if "`var'"=="scale" {
		replace col = " $ iva_{t} $" if [_n]==`i'
	}
	else if "`var'"=="assetgrowth" {
		replace col = " $ ag_{t} $" if [_n]==`i'
	}	
	else {
		replace col = " $ `var'_{t} $" if [_n]==`i'
	}
	local i = `i' + 1
}
 
* Keep only state variables
replace col = "\addlinespace " + col if [_n]>=2
order col
foreach var of varlist B* {
	qui tostring `var', force replace format("%5.3f")
} 

local c = ""
local columns = "" 
foreach var in `state_vars' {
	if "`var'"=="scale" {
		local columns = "`columns' & $ iva_{t} $"
	}
	else if "`var'"=="assetgrowth" {
		local columns = "`columns' & $ ag_{t} $"
	}	
	else{
		local columns = "`columns' & $ `var'_{t} $"
	}
	local c = "`c' c"
}	
		 
listtab * using "Results/Tables/02_15_`vweight'_`end'.tex", ///
	rstyle(tabular) replace ///
	head("\begin{tabular}{l`c'}" ///
	"\midrule" ///
	"`columns' \\" ///
	"\midrule") ///
	foot("\midrule" "\end{tabular}")	
	
****************************************************************************
* 02.16  Decomposition of return news   (Table 02_16_`vweight')
**************************************************************************** 

mat B = (B \ (1,0,0,0,0,0,0,0,0,0,0,0,0)) /* r_L1 */
mat B = (B \ (0,1,0,0,0,0,0,0,0,0,0,0,0)) /* roe_L1 */
mat B = (B \ (0,0,1,0,0,0,0,0,0,0,0,0,0)) /* iva_L1 */
mat B = (B \ (0,0,0,0,1,0,0,0,0,0,0,0,0)) /* ag_L1 */
mat B = (B \ (0,0,0,0,0,1,0,0,0,0,0,0,0)) /* r_L2 */
mat B = (B \ (0,0,0,0,0,0,0,0,1,0,0,0,0)) /* ag_L2 */
mat B = (B \ (0,0,0,0,0,0,0,0,0,1,0,0,0)) /* r_L3 */
mat B = (B \ (0,0,0,0,0,0,0,0,0,0,1,0,0)) /* ag_L3 */

preserve

clear
svmat B
	* Save the transition matrix in csv
	export excel "Data/created/extended_VAR.xlsx", replace

restore

local num_vars = 13
matrix Sigma_new = J(`num_vars',`num_vars',0)
local N = 5
forvalues i = 1/`N' {
	forvalues j = 1/`N' {
		mat Sigma_new[`i',`j'] = Sigma[`i',`j']
	}
} 
mat Sigma = Sigma_new

* Define elementary vectors  
forvalues i = 1/3 {
	matrix e`i' = J(`num_vars',1,0) 
	matrix e`i'[`i',1] = 1 
} 

* Generate lambdas 
matrix lambda =	(e1'*`rho'*B*inv(I(`num_vars')-`rho'*B))'
matrix lambda2 = (e2'*inv(I(`num_vars')-`rho'*B))'						// For computing roe news explicitly
matrix lambda3 = (e3'*inv(I(`num_vars')-`rho'*B))'						// For computing scale news explicitly
matrix lambda2r = (e1 + lambda - lambda3)								// For computing roe news as a residual term
matrix lambda3r = (e1 + lambda - lambda2)								// For computing scale news as a residual term

* Generate news variance and covariance
* Back out everything explicitly
* Variances 
matrix var_N_DR = lambda'*Sigma*lambda
matrix var_N_roe = lambda2'*Sigma*lambda2
matrix var_N_scale = lambda3'*Sigma*lambda3

* Covariances	
matrix cov_N_DR_roe = lambda'*Sigma*lambda2
matrix cov_N_DR_scale = lambda'*Sigma*lambda3
matrix cov_N_roe_scale = lambda2'*Sigma*lambda3	


* Variance-covariance matrix	
local var var_N_DR var_N_roe var_N_scale cov_N_DR_roe cov_N_DR_scale cov_N_roe_scale
foreach mat in `var' {
	local `mat' = `mat'[1,1]
}
matrix totalvar = var_N_DR + var_N_roe + var_N_scale - 2*(cov_N_DR_roe + cov_N_DR_scale - cov_N_roe_scale)
local totalvar = totalvar[1,1]
local svarDR = `var_N_DR' / `totalvar'
local svarroe = `var_N_roe' / `totalvar'
local svarscale = `var_N_scale' / `totalvar'
local svarDRroe = -2*`cov_N_DR_roe' / `totalvar'
local svarDRscale = -2*`cov_N_DR_scale' / `totalvar'
local svarroescale = 2*`cov_N_roe_scale' / `totalvar'			
local var var_N_DR var_N_roe var_N_scale
foreach mat in `var' {
	local `mat' = sqrt(`mat'[1,1])
}
local cov_N_DR_roe = `cov_N_DR_roe'/(`var_N_DR'*`var_N_roe')
local cov_N_DR_scale = `cov_N_DR_scale'/(`var_N_DR'*`var_N_scale')
local cov_N_roe_scale = `cov_N_roe_scale'/(`var_N_roe'*`var_N_scale')	

mat VCM = (`var_N_DR', 0, 0, 0, `svarDR', 0, 0 \ ///
		   `cov_N_DR_roe', `var_N_roe', 0,0, `svarDRroe', `svarroe', 0 \ ///
		   `cov_N_DR_scale', `cov_N_roe_scale', `var_N_scale', 0, `svarDRscale', `svarroescale', `svarscale')

* Export to Latex 
clear
svmat VCM
g col = ""
g row = 0
local i = 1
foreach var in "N_{DR}" "N_{roe}" "N_{iva}" {
		replace col = "$ `var' $" if [_n]==`i'
		replace row = `i'-0.5 if [_n]==`i'
		local i = `i' + 1
		}
replace col = "\addlinespace " + col if [_n]>=2
order col 

*			append using "Output/news_decomp_stderr_C.dta"	

forval i= 1/7 {
	qui tostring VCM`i', force replace format("%5.3f")
	qui replace VCM`i' = "" if VCM`i' == "0.000" 
}

sort row
drop row

* Output to latex
listtab * using "Results/Tables/02_16_`vweight'_`end'.tex", ///
	rstyle(tabular) replace ///
	head("\begin{tabular}{lccccccc}" ///
	"\midrule" ///
	"& \multicolumn{3}{c}{$\sigma$ (diag), $\rho$ (off-diag)} & & \multicolumn{3}{c}{Contribution to $\sigma^2_{N_r}$} \\" ///
	"\cmidrule(lr){2-4} \cmidrule(lr){6-8}" /// 
		" & $ N_{DR} $ & $ N_{roe} $ & $ N_{iva} $ & & $ -N_{DR} $ & $ N_{roe} $ & $ N_{iva} $\\" ///
	"\midrule") ///
	foot("\midrule" "\end{tabular}")
 
****************************************************************************
* 02.17  Decompose the cross-sectional var(mb)  (Table 02_17_`vweight')
**************************************************************************** 	 
use "`resid'", clear
local lags2=`lags'-1
local nonr_vars = "roe scale bm assetgrowth"		
forval i = 1/`lags2' {
	local nonr_vars = "`nonr_vars' r_L`i'"
}
foreach var in roe scale assetgrowth {
	forval i = 1/`lags2' {
		if "`var'"=="assetgrowth"{
			local nonr_vars = "`nonr_vars' `var'_L`i'" 
		}
		else if inlist("`var'","roe","scale","sgrowth","gamma") & `i'==1 {
		local nonr_vars = "`nonr_vars' `var'_L`i'" 
		}
	}
}
	matrix expected = inv(I(`num_vars')-`rho'*B) * B
	mat list expected
	matrix expected_r = expected[1, 1..`num_vars']
	matrix expected_roe = expected[2, 1..`num_vars'] 
	matrix expected_scale = expected[3, 1..`num_vars']
	
	foreach var in r roe scale {
		svmat expected_`var'
		forvalues i = 1/`num_vars' {
			replace expected_`var'`i' = expected_`var'`i'[_n-1] if [_n]>1
		}
		gen exp_`var' = r * expected_`var'1
		local i = 2
		foreach var2 in `nonr_vars' {
			replace exp_`var' = exp_`var' + `var2' * expected_`var'`i'
			local i = `i' + 1
		}
		drop expected_`var'*
	} 

	* Cross-sectional decomposition
	local i = 1
	replace exp_r = -1 * exp_r
	foreach var in r roe scale {
		label var mb "$ mb_t $"
		reghdfe exp_`var' mb [aweight = weight], nocon absorb(year) cl(year firm)  
			* Save results	  
			eststo C`i'
			local i = `i' + 1	
	}
	// Tentative table in latex
	esttab C1 C2 C3 using ///
		"Results/Tables/02_17_`vweight'_`end'.tex", ///
		 b(3) se /*t(1)*/  margin /*booktabs*/ nomtitles ///
		 nostar label ///
		 mgroups("$ -\sum_{j=1}^{\infty}\rho^j E_t\left[r_{t+j}\right] $" ///
			"$ \sum_{j=1}^{\infty} E_t\left[roe_{t+j}\right] $" ///
			"$ \sum_{j=1}^{\infty} E_t\left[iva_{t+j}\right] $", pattern(1 1 1)                   ///
		 prefix(\multicolumn{@span}{c}{) suffix(})   ///
		 span erepeat(\cmidrule(lr){@span}))    /// 
		mlabels(none) collabels(none) stats(N, fmt(%5.0fc) ///
		labels("Observations")) /// 
		nonotes replace constant substitute(\_ _) 		 
		
	eststo clear

****************************************************************************
* 02.18  Predictability  (Figure 4)
**************************************************************************** 
* Choose window length in years
local window = 5
local start = 1965																 
local end = 2022
local rho = 0.96
********************************************************************************
* 02.18.01  Load data
******************************************************************************** 
* Load annual dataset
u "Data/Created/Annual Data 2021 with Bins", clear 
rename scale iva
rename scale_L1 iva_L1
gen cf = roe + iva


drop if year<`start'

********************************************************************************
* 02.18.02  Cross-sectional univariate forecasting
******************************************************************************** 

* Set up Window

qui sum year, d
local begy = `r(min)' + `window'
local endy = `r(max)'

* List of LHS variables

local list = "r roe iva cf mb"

* Matrix to store coefficients, SE's

foreach var in `list'{
	mat B`var' = J(1,1,.)
	mat SE`var' = J(1,1,.)
	mat VAR_err`var' = J(1,1,.)
	mat cov_mb_`var'_impl = J(1,1,.)
}

mat phiCS = J(1,1,.)
mat NoObs = J(1,1,.)
mat SDmb_L1 = J(1,1,.)
mat VAR_err_mb_L1_impl = J(1,1,.)

* Demean the mb_L1

qui egen aux = mean(mb_L1), by(year)
qui g demean_mb_L1 = aux - mb_L1
drop aux

* Univariate regressions forecasting 1-year return

forval i = `begy'/`endy'{
	foreach var of varlist `list'{
	
		* Cross-sectional retression
		qui reghdfe `var' mb_L1 if year < `i' & year >= `i' - `window', ///
			a(year) cl(year gvkey) resid
	
		* Store coefficients and standard errors
		mat B`var' = (B`var' \ _b[mb_L1])
		mat SE`var' = (SE`var' \ _se[mb_L1])
		
		* Residual
		qui predict err_`var', resid

	}

	* Cross-sectional autocorrelations
	qui reghdfe mb mb_L1 if year < `i' & year >= `i' - `window', a(year) cl(year gvkey)
	mat phiCS = (phiCS \ _b[mb_L1])
	

	
	* Store the number of observations
	mat NoObs = (NoObs \ e(N))
	
	* Standard deviation of mb_L1 after demeaning 
	qui sum demean_mb_L1 if year < `i' & year >= `i' - `window', d
	mat SDmb_L1 = (SDmb_L1 \ r(sd))
	
	* Implied residual of mb autocorrelation
	g err_mb_impl = (err_r - err_cf)/`rho' 
	qui sum err_mb_impl, d
	mat VAR_err_mb_L1_impl = (VAR_err_mb_L1_impl \ r(sd)^2)
	
	* Variance-covariance matrix with the implied estimates of e_{mb}
	foreach var of varlist `list'{
		qui corr err_mb_impl err_`var', cov
		mat VAR_err`var' = (VAR_err`var' \ r(Var_2))
		mat cov_mb_`var'_impl = (cov_mb_`var'_impl \ r(cov_12))
		}

	drop err_*
}
	
* Collect stored estimates 

clear
foreach var in `list' {
	svmat B`var'
	svmat SE`var'
	svmat VAR_err`var'
	svmat cov_mb_`var'_impl
}
svmat phiCS
svmat NoObs
svmat SDmb_L1
svmat VAR_err_mb_L1_impl
drop if _n == 1

* Rename variables 

foreach var in `list' {
	ren B`var'1 B`var'
	ren SE`var'1 SE`var'
	ren VAR_err`var'1 VAR_err`var'
	ren cov_mb_`var'_impl1 cov_mb_`var'_impl
}

ren phiCS1 phiCS
ren NoObs1 NoObs
ren SDmb_L11 SDmb_L1
ren VAR_err_mb_L1_impl1 VAR_err_mb_L1_impl
replace NoObs = NoObs - 1 	/* Correct for degree-of-freedom */

* Implied autocorrelation estimates

g phiCS_impl = (1 + Br - Bcf)/`rho'

* Compute t-stats and long-run coefficients
foreach var in  `list' {
	gen VAR_err`var'2 = SE`var'^2*NoObs
	gen b_LRCS_`var' = B`var'/(1-(`rho' * phiCS))
	gen t_`var' = B`var' / SE`var'
	gen SE_LRCS_`var' = 1/(1 - `rho'*phiCS_impl)* ///
		(VAR_err`var'2 + 2*`rho'*b_LRCS_`var'*cov_mb_`var'_impl ///
			+ `rho'^2*b_LRCS_`var'^2*VAR_err_mb_L1_impl)^0.5 ///
		/(NoObs^0.5*SDmb_L1)
	gen SE_SRCS_`var' = VAR_err`var'^0.5 / NoObs^0.5	
	}

gen checkCS = b_LRCS_cf - b_LRCS_r

replace Br = -Br
replace b_LRCS_r = -b_LRCS_r

* Set year

g year = `begy'
replace year = year[_n-1] + 1 if _n > 1
tsset year

* Set confidence interval
foreach var in `list' {
	g up_B`var' = B`var' + 1.96*SE`var'
	g down_B`var' = B`var' - 1.96*SE`var'
	g up_B`var'2 = B`var' + 1.96*(VAR_err`var'/NoObs)^0.5
	g down_B`var'2 = B`var' - 1.96*(VAR_err`var'/NoObs)^0.5
	g up_LRCS_`var' = b_LRCS_`var' + 1.96*SE_LRCS_`var'
	g down_LRCS_`var' = b_LRCS_`var' - 1.96*SE_LRCS_`var'	
	g up_LRCS_`var'2 = b_LRCS_`var' + 1.96*SE_LRCS_`var'*SE`var'/(VAR_err`var'/NoObs)^0.5
	g down_LRCS_`var'2 = b_LRCS_`var' - 1.96*SE_LRCS_`var'*SE`var'/(VAR_err`var'/NoObs)^0.5	
}
	
* Label variables

label var b_LRCS_r  		"LR r"
label var b_LRCS_mb  		"LR mb"
label var b_LRCS_roe  		"LR roe"
label var b_LRCS_iva  		"LR iva"
label var b_LRCS_cf  		"LR cf"

label var Br  		"r"
label var Bmb  		"mb"
label var Bcf  		"cf"
label var Broe 		"roe"
label var Biva 		"iva"

* Plot LR with clustered delta-method CIs
twoway  (line Bcf year, lcolor(maroon) lpattern(dash) yaxis(1)) ///
		(line Br year, lcolor(navy) lpattern(dash) yaxis(1)) ///
		(rarea up_LRCS_cf down_LRCS_cf year, fcolor(red%20) lwidth(none) yaxis(2)) ///
		(line b_LRCS_cf year, lcolor(red) yaxis(2)) ///
		(rarea up_LRCS_r down_LRCS_r year, fcolor(gs4%20) lwidth(none) yaxis(2)) ///
		(line b_LRCS_r year, lcolor(gs4) yaxis(2)), ///
		ttitle("") graphregion(color(white) lcolor(white)) ylabel(-0.2(0.1)0.3, axis(1)) yscale(range(-0.2 0.3) axis(1)) ylabel(-1(1)2, axis(2)) yscale(range(-1.333 2) axis(2)) ///
		legend(order(1 2 4 6) label(1 "SR: Cash flows (left)") label(2 "SR: Discount rates (left)") label(4 "LR: Cash flows (right)") label(6 "LR: Discount rates (right)"))

gr export "Results/Figures/02_18.pdf", replace as(pdf)		
